#include "StdAfx.h"
#include "u_i.h"

u_i::u_i(void)
{
}
u_i::u_i(double Rn,double Ic,double C, double w):rscjprops(Rn,Ic,C,w),utobj(Rn,Ic,C,0,0,0,w,0.01){
	this->Iac=0;
	this->iac=0;
	this->Imin=0;
	this->imin=0;
	this->Imax=5.0/Ic;
	this->imax=5;
	this->deltaI=0.01*Ic;
	this->deltai=this->deltaI/this->Ic;
	this->h=0.01;
	this->calcbeta();
	std::cout<<this->betac<<std::endl;
	this->calcomegac();
	this->calcomega();
	this->filename="data";
	this->constI=true;
	this->Iasy=0;
	this->iasy=0;
	this->Iex=0;
	this->iex=0;
}
u_i::u_i(double betac , double omega):rscjprops(betac,omega),utobj(betac,0,0,0,omega,0.01){
	this->Iac=0;
	this->iac=0;
	this->imax=5;
	this->Imax=-1;
	this->imin=0;
	this->Imin=-1;
	this->deltai=0.01;
	this->deltaI=-1;
	this->h=0.01;
	this->betac=betac;
	this->omega=omega;
	this->filename="data";
	this->constI=true;
	this->Iasy=-1;
	this->iasy=0;
	this->Iex=-1;
	this->iex=0;
}

void u_i::dec_dI(){
	this->I.reserve(this->I.size()*2);
	 std::vector<long double>::iterator it;
	 it=this->I.begin();
	for(int i=0;i<this->I.size()-1;i+=2){
		I.insert(it+i+1,(I.at(i)+I.at(i+1))/2);
	}
}

void u_i::inv_dec_dI(){
	std::vector<long double> Itemp;
	std::vector<long double> Utemp;
	Itemp.reserve(this->I.size()/2);
	Utemp.reserve(this->U.size()/2);
	for(int i=0;i<this->I.size();i+=2){
		Itemp.push_back(this->I.at(i));
		Utemp.push_back(this->U.at(i));
	}
	this->I=Itemp;
	this->U=Utemp;
}

void u_i::work(){
		
	if(this->constI){
		I.clear();
		U.clear();
		I.resize(0);
		U.resize(0);

		int n=2*((imax-imin)/deltai)+4;
		I.reserve(n);
		U.reserve(n);
		I.push_back(0);

		utobj.setiac(iac);
		U.push_back(calcu(0,0));
	int i=0;
		while(I[i]<imax){
			++i;
			I.push_back(I[i-1]+deltai);
			U.push_back(calcu(I[i],U[i-1]));
		}
		while(I[i]>imin){
			++i;
			I.push_back(I[i-1]-deltai);
			U.push_back(calcu(I[i],U[i-1]));
		}
		while(I[i]<0){
			++i;
			I.push_back(I[i-1]+deltai);
			U.push_back(calcu(I[i],U[i-1]));
		}
	}else{
		/*if(this->T!=0){
		this->dec_dI();
		}*/
		U.reserve(I.size());
		this->imax=*max_element(I.begin(),I.end());
		this->imin=*min_element(I.begin(),I.end());
		utobj.setiac(iac);
		U.push_back((I.at(0)));
		for(int i=1;i<I.size();++i){
			U.push_back(calcu(I[i],U[i-1]));
		}
	}
	
	
		
	
	//write();
	if(!normalized){
		for(int i=0;i<I.size();++i){
			I.at(i)*=this->Ic;
			U.at(i)*=(this->Ic*this->Rn);
		}
	
	}
	if((this->T>0)&&(!normalized)){
		this->temp();

		/*if(!(this->constI)){
			this->inv_dec_dI();
		}*/
	}


}

long double u_i::calcu(long double i,long double u){
	long double U1=0;
	double itemp;
	itemp=i-this->iasy;
	if(itemp>0){
		itemp=itemp-this->iex;
	}else if(itemp<0){
		itemp=itemp+this->iex;
	}
		utobj.seth(h);
		utobj.setbetac(this->betac);
		utobj.setomega(this->omega);
		utobj.setiac(this->iac);
		utobj.setidc(itemp);
		utobj.setu0(u);
		utobj.setsmax(2000);
	utobj.work();
	U1=utobj.getu();
	return U1;
}





void u_i::write(){
	
	std::fstream datei;
	datei.open(filename, std::ios::out);
	if(normalized){
		datei<<"#Revision="<<REVISION<<std::endl;
		datei<<"#Betac="<<betac<<std::endl;
		datei<<"#Omega="<<omega<<std::endl;
		datei<<"#Iac="<<iac<<std::endl;
		datei<<"#Iasy="<<iasy<<std::endl;
		datei<<"#Iex="<<iex<<std::endl;
	}else{
		datei<<"#Revision="<<REVISION<<std::endl;
		datei<<"#Betac="<<betac<<std::endl;
		datei<<"#Omega="<<omega<<std::endl;
		datei<<"#Omegac="<<wc<<std::endl;
		datei<<"#Rn="<<Rn<<"Ohm"<<std::endl;
		datei<<"#C="<<C<<"nF"<<std::endl;
		datei<<"#Ic="<<Ic<<"mA"<<std::endl;
		datei<<"#Iac="<<Iac<<"mA"<<std::endl;
		datei<<"#Iasy="<<Iasy<<"mA"<<std::endl;
		datei<<"#Iex="<<Iex<<"mA"<<std::endl;
	}
	int n=I.size();
	for(int i=0;i<n;i++){
		datei<<I[i]<<" "<<U[i]<<std::endl;
	}
	datei.close();

	
	/*datei.open("gaus.txt", std::ios::out);
	for(int i=0;i<this->irnd.size();i++){
		datei<<this->irnd.at(i)<<std::endl;
	}
	datei.close();*/

}
u_i::~u_i(void)
{
}

void u_i::setimax(double imax){
	this->normalized=true;
	this->imax=imax;
	this->Imax=-1;
}
void u_i::setImax(double Imax){
	this->normalized=false;
	this->Imax=Imax;
	if(Ic>0){
		this->imax=Imax/Ic;
	}
}
void u_i::setimin(double imin){
	this->normalized=true;
	this->imin=imin;
	this->Imin=-1;
}
void u_i::setImin(double Imin){
	this->normalized=false;
	this->Imin=Imin;
	if(Ic>0){
		this->imin=Imin/Ic;
	}
}
void u_i::setiasy(double iasy){
	this->normalized=true;
	this->iasy=iasy;
	this->Iasy=-1;
}
void u_i::setIasy(double Iasy){
	this->normalized=false;
	this->Iasy=Iasy;
	if(Ic>0){
		this->iasy=Iasy/Ic;
	}
}
void u_i::setiex(double iex){
	this->normalized=true;
	this->iex=iex;
	this->Iex=-1;
}
void u_i::setIex(double Iex){
	this->normalized=false;
	this->Iex=Iex;
	if(Ic>0){
		this->iex=Iex/Ic;
	}
}

void u_i::setDeltai(double Deltai){
	this->normalized=true;
	this->constI=true;
	this->deltai=Deltai;
	this->deltaI=-1;
}
void u_i::setDeltai(std::vector<double> Deltai){
	this->normalized=true;
	this->constI=false;
	this->I.clear();
	for(int i=0;i<Deltai.size();++i){
		this->I.push_back(Deltai.at(i));
	}
}

void u_i::setDeltaI(std::vector<double> DeltaI){
	this->normalized=false;
	this->constI=false;
	this->I.clear();
	for(int i=0;i<DeltaI.size();++i){
		this->I.push_back(DeltaI.at(i)/this->Ic);
	}
}

void u_i::setDeltaI(double DeltaI){
	this->normalized=false;
	this->constI=true;
	this->deltaI=DeltaI;
	if(Ic>0){
		this->deltai=DeltaI/Ic;
	}
}
void u_i::setsmax(double smax){
	this->smax=smax;
	this->utobj.setsmax(smax);
}
void u_i::setfilename(std::string name){
this->filename=name;

}

double u_i::getimin(){
	return imin;
}
double u_i::getImin(){
	return Imin;
}
double u_i::getiasy(){
	return iasy;
}
double u_i::getIasy(){
	return Iasy;
}
double u_i::getiex(){
	return iex;
}
double u_i::getIex(){
	return Iex;
}
double u_i::getimax(){
	return imax;
}
double u_i::getImax(){
	return Imax;
}
double u_i::getDeltai(){
	return deltai;
}
double u_i::getDeltaI(){
	return Imax;
}
double u_i::getsmax(){
	
	return utobj.getsmax();
}


std::string u_i::getfilename(){
	return this->filename;
}

double getmax(double *x,int n){
	double temp=x[0];
	for(int i=1;i<n;++i){
		if(temp<x[i]){
			temp=x[i];
	}
		
}
	return temp;
}


double getmin(double *x,int n){
	double temp=x[0];
	for(int i=1;i<n;++i){
		if(temp>x[i]){
			temp=x[i];
	}
		
}
	return temp;
}

double u_i::interpol_lin(std::vector<long double> xdata,std::vector<long double> ydata, double x,int n1){
	int i=n1;
	double *xd=0;
	double *yd=0;
	double temp;
	int sig;
	int n=xdata.size()/2;
	int n2=xdata.size()-n;

	if(xdata.at(0)==xdata.back()){
		
		if(i<n){

			if(xdata.at(0)>xdata.at(1)){
				sig=-1;
			}else{
				sig=1;
			}
			xd = new double [n];
			yd = new double [n];
			for(int i=0;i<n;++i){
				xd[i]=sig*xdata.at(i);
				yd[i]=ydata.at(i);
			}
		}else{

			if(xdata.at(n)>xdata.at(n+1)){
				sig=-1;
			}else{
				sig=1;
			}
			xd = new double [n2];
			yd = new double [n2];
			for(int i=n;i<xdata.size();++i){
				xd[i-n]=sig*xdata.at(i);
				yd[i-n]=ydata.at(i);
			}

		}
	}else{
		n=xdata.size();
		if(xdata.at(0)>xdata.at(1)){
				sig=-1;
			}else{
				sig=1;
			}
			xd = new double [n];
			yd = new double [n];
			for(int i=0;i<n;++i){
				xd[i]=sig*xdata.at(i);
				yd[i]=ydata.at(i);
			}

	}	
	
	if((sig*x<getmin(xd,n))||(sig*x>getmax(xd,n))){
		return ydata.at(i);
	}
	

	
	gsl_interp_accel *acc = gsl_interp_accel_alloc ();
	/*gsl_spline *spline = gsl_spline_alloc (gsl_interp_linear, n);
	gsl_spline_init (spline, xd, yd, n);
	temp=gsl_spline_eval (spline, sig*x, acc);*/
	//gsl_spline_free (spline);
	//gsl_interp_accel_free (acc);


	gsl_interp *interp= gsl_interp_alloc (gsl_interp_linear, n);
	gsl_interp_init (interp, xd, yd, n);
	temp=gsl_interp_eval (interp, xd, yd, sig*x, acc);
	gsl_interp_free (interp);
	gsl_interp_accel_free (acc);
	//delete spline;
	//delete acc;
	delete xd;
	delete yd;
	return temp;
}

void u_i::temp(){

	//this->acc = gsl_interp_accel_alloc ();
	//this->spline = gsl_spline_alloc (gsl_interp_linear, this->I.size());
	double sigma=sqrt((2.0*1.3806488*this->T)/(PI*pow(this->Rn,2)*this->C))*1e-4;
	std::cout<<"sigma="<<sigma<<std::endl;
	std::vector<long double> Utemp;
	Utemp.resize(this->U.size());
	
	/*const gsl_rng_type * type;

	gsl_rng_env_setup();
	type = gsl_rng_default;
	gsl_rng * r = gsl_rng_alloc (gsl_rng_ranlxd2);
	gsl_rng_set(r,0);*/


	for(int i=0;i<Utemp.size();++i){
		Utemp.at(i)=0;
	}

	int count=500;

	/*for(int i=0;i<count;++i){
		this->irnd.push_back(gsl_ran_gaussian (r, sigma));
	}*/

	this->irnd.clear();
	double deltaItemp=0.1;
	double A=0.5*count*pow(sigma,2);
	long double exponent;
	bool temp=true;
	int ntemp=0;
	while(temp){
		for(int n=-(5.0/deltaItemp);n<(5.0/deltaItemp);++n){
			exponent=-0.5*(pow(((n*deltaItemp)),2.0));
			ntemp=ceil(20.0*exp(exponent)-0.5);
			
			if(ntemp>0){
			for(int k=0;k<ntemp;++k){
				this->irnd.push_back(n*deltaItemp*sigma);
			}
		}
		}
			temp=false;
	}
	count=this->irnd.size();
	for(int i=0;i<this->I.size();++i){
		
		
		for(int k=0;k<count;++k){
			
			Utemp.at(i)+=interpol_lin(this->I,this->U,this->I.at(i)+this->irnd.at(k),i);
		}
		Utemp.at(i)=Utemp.at(i)/count;
	}
	
	this->U=Utemp;
	Utemp.~vector();
	//gsl_rng_free (r);
	//gsl_spline_free (this->spline);
	//gsl_interp_accel_free (this->acc);
}


